# Reproduction files - Study 1 descriptives, opinion formation

# Deepening Bridging and Moving Minds in Stressful Times







# Session Information -----------------------------------------------------




#R version 4.2.2 (2022-10-31 ucrt)
#Platform: x86_64-w64-mingw32/x64 (64-bit)
#Running under: Windows 10 x64 (build 19045)
#
#Matrix products: default
#
#locale:
#  [1] LC_COLLATE=German_Germany.utf8  LC_CTYPE=German_Germany.utf8    LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C                   
#[5] LC_TIME=German_Germany.utf8    
#
#attached base packages:
#  [1] stats     graphics  grDevices utils     datasets  methods   base     
#
#other attached packages:
#  [1] car_3.1-1       carData_3.0-5   patchwork_1.3.0 scales_1.3.0    ggeffects_1.1.5 stargazer_5.2.3 nnet_7.3-18     sjPlot_2.8.12  
#[9] psych_2.2.9     forcats_1.0.0   stringr_1.5.0   dplyr_1.1.4     purrr_1.0.1     readr_2.1.4     tidyr_1.3.0     tibble_3.2.1   
#[17] ggplot2_3.5.1   tidyverse_1.3.2 haven_2.5.3     sjmisc_2.8.9   
#
#loaded via a namespace (and not attached):
#  [1] nlme_3.1-160        fs_1.6.1            lubridate_1.9.2     insight_0.19.0      httr_1.4.4          tools_4.2.2         backports_1.4.1    
#[8] R6_2.5.1            sjlabelled_1.2.0    DBI_1.1.3           colorspace_2.1-0    withr_2.5.0         tidyselect_1.2.0    Exact_3.2          
#[15] mnormt_2.1.1        emmeans_1.8.4-1     compiler_4.2.2      performance_0.10.2  cli_3.6.0           rvest_1.0.3         expm_0.999-7       
#[22] xml2_1.3.3          sandwich_3.0-2      bayestestR_0.13.0   mvtnorm_1.1-3       proxy_0.4-27        minqa_1.2.5         pkgconfig_2.0.3    
#[29] lme4_1.1-31         dbplyr_2.3.0        rlang_1.1.4         readxl_1.4.2        rstudioapi_0.14     farver_2.1.1        generics_0.1.3     
#[36] zoo_1.8-11          jsonlite_1.8.4      googlesheets4_1.0.1 magrittr_2.0.3      Matrix_1.5-1        Rcpp_1.0.10         DescTools_0.99.50  
#[43] munsell_0.5.0       abind_1.4-5         lifecycle_1.0.3     stringi_1.7.12      multcomp_1.4-22     MASS_7.3-58.1       rootSolve_1.8.2.3  
#[50] grid_4.2.2          parallel_4.2.2      crayon_1.5.2        lmom_2.9            lattice_0.20-45     splines_4.2.2       sjstats_0.18.2     
#[57] hms_1.1.2           knitr_1.42          pillar_1.10.1       ggpubr_0.6.0        boot_1.3-28         gld_2.6.6           estimability_1.4.1 
#[64] ggsignif_0.6.4      codetools_0.2-18    reprex_2.0.2        glue_1.6.2          data.table_1.14.6   modelr_0.1.10       nloptr_2.0.3       
#[71] vctrs_0.6.5         tzdb_0.3.0          cellranger_1.1.0    gtable_0.3.1        assertthat_0.2.1    datawizard_0.6.5    xfun_0.37          
#[78] xtable_1.8-4        broom_1.0.3         e1071_1.7-13        coda_0.19-4         rstatix_0.7.2       survival_3.4-0      class_7.3-20       
#[85] googledrive_2.0.0   gargle_1.3.0        timechange_0.2.0    TH.data_1.1-1       ellipsis_0.3.2 


# Calculations ------------------------------------------------------------







# load necessary packages
library(sjmisc)
library(haven)
library(tidyverse)
library(psych)
library(sjPlot)
library(nnet)
library(stargazer)
library(ggeffects)
library(scales)
library(patchwork)
library(car)


# 1. Data preperation, descriptives ---------------------------------------

# read dataset
data_DE <- readRDS("data_de_wght.RDS")

# number of cases
nrow(data_DE) 
frq(data_DE$finished_post) # 2132 completes in germany

# descriptive information pre and post survey position
describe(data_DE$position_pre_cont) # pre: 6.71
describe(data_DE$position_post_cont) # post: 6.7

# create new variable: pro health vs. pro liberties
data_DE$pos_pre_group <- car::recode(data_DE$position_pre_cont,
                                     recodes = "1:5 = 0; 6:10 = 1", as.factor = T)
frq(data_DE$pos_pre_group)
# 32.09 pro freedom; 67.91 pro health

# new variable strong prior opinion
data_DE$strong_prior <- car::recode(data_DE$position_pre_cont,
                                    recodes = "c(1,2,9,10) = 1; 3:8 = 0", as.factor = T)
frq(data_DE$strong_prior)
# 37.32 % strong prior; 62.68 % no strong prior


# diversity of the sample

# age
describe(data_DE$age) # mean = 46.08, min = 18, max = 89

# gender
frq(data_DE$gender) # 48.26 female, 51.50 male, 0.23 diverse

# recode education (low, mid, high)
data_DE$education_rec <- car::recode(data_DE$education,
                                     recodes = "1:3=1; 4=2; 5:7= 3", as.factor = T)
frq(data_DE$education_rec) # 11.07 low; 30.82 medium; 58.11 high
frq(data_DE$education)


# Plot pre position dist study 1 & 2

p1 <- ggplot(data_DE, aes(position_pre_cont))+
  geom_histogram(bins = 10, fill = c(rep("red",2), rep("grey",6), rep("red",2)))+
  theme_bw()+
  ggtitle("Study 1: Germany (March 2021)")+
  xlab("Position Pre Survey")+
  ylab("Count")+
  geom_vline(xintercept = 5.5, alpha = 0.5)+
  geom_vline(aes(xintercept=mean(position_pre_cont, na.rm =T)), col = "red", linetype = "dashed")+
  annotate("text", x = mean(data_DE$position_pre_cont, na.rm =T)-0.1, y = attr(DescTools::Mode(data_DE$position_pre_cont, na.rm =T), "freq"), label = paste0("Mean = ", round(mean(data_DE$position_pre_cont, na.rm =T),1)), hjust = 1)+
  theme(legend.position = "none")+
  scale_x_continuous(breaks = c(1:10))+
  ylim(0,550)

# read AT dataset
data_AT <- readRDS("data_at_wght.RDS")

p2 <- ggplot(data_AT, aes(position_pre_cont))+
  geom_histogram(bins = 10, fill = c(rep("red",2), rep("grey",6), rep("red",2)))+
  theme_bw()+
  ggtitle("Study 2: Austria (December 2021)")+
  xlab("Position Pre Survey")+
  ylab("Count")+
  geom_vline(xintercept = 5.5, alpha = 0.5)+
  geom_vline(aes(xintercept=mean(position_pre_cont, na.rm =T)), col = "red", linetype = "dashed")+
  annotate("text", x = mean(data_AT$position_pre_cont, na.rm =T)-0.1, y = attr(DescTools::Mode(data_AT$position_pre_cont, na.rm =T), "freq"), label = paste0("Mean = ", round(mean(data_AT$position_pre_cont, na.rm =T),1)), hjust = 1)+
  theme(legend.position = "none")+
  scale_x_continuous(breaks = c(1:10))+
  ylim(0,550)





p1 <- ggplot(data_DE, aes(position_pre_cont))+
  geom_histogram(bins = 10, fill = "darkgrey")+
  theme_bw()+
  ggtitle("Study 1: Germany (March 2021)")+
  xlab("Position Pre Survey")+
  ylab("Count")+
  geom_vline(aes(xintercept=mean(position_pre_cont, na.rm =T)), col = "red", linetype = "dashed")+
  annotate("text", x = mean(data_DE$position_pre_cont, na.rm =T)-0.1, y = 500, label = paste0("Mean = ", round(mean(data_DE$position_pre_cont, na.rm =T),1)), hjust = 1)+
  theme(legend.position = "none")+
  scale_x_continuous(breaks = c(1:10))+
  ylim(0,760)

p2 <- ggplot(data_AT, aes(position_pre_cont))+
  geom_histogram(bins = 10, fill = "darkgrey")+
  theme_bw()+
  ggtitle("Study 2: Austria (December 2021)")+
  xlab("Position Pre Survey")+
  ylab("Count")+
  geom_vline(aes(xintercept=mean(position_pre_cont, na.rm =T)), col = "red", linetype = "dashed")+
  annotate("text", x = mean(data_AT$position_pre_cont, na.rm =T)-0.1, y = 500, label = paste0("Mean = ", round(mean(data_AT$position_pre_cont, na.rm =T),1)), hjust = 1)+
  theme(legend.position = "none")+
  scale_x_continuous(breaks = c(1:10))+
  ylim(0,760)


p_comb <- p1+p2

ggsave("position_pre.png", plot = p_comb, width = 10, height = 6)




# 2. Opinion change -------------------------------------------------------

# descriptive information pre and post survey opinion
describe(data_DE$position_pre_cont) # mean pre: 6.71
describe(data_DE$position_post_cont) # mean post: 6.70

# aggregate opinion change
data_DE <- data_DE %>%
  mutate(abs_position_change = abs(data_DE$position_post_cont - data_DE$position_pre_cont),
         position_change1 = ifelse(abs_position_change <= 0, "no change", "change"),
         position_change2 = ifelse(abs_position_change <=1, "no change", "change"))

# average individual opinion shift
mean(data_DE$abs_position_change, na.rm = T) # avg. change 0.75 scale points




# calculate abs diff
data_DE$abs_change = abs(data_DE$position_post_cont - data_DE$position_pre_cont)

# mean
mean_abs_change <- mean(data_DE$abs_change, na.rm = TRUE)

# sd
sd_abs_change <- sd(data_DE$abs_change, na.rm = TRUE)

# Cohen's d
cohen_d_abs <- mean_abs_change / sd_abs_change
cohen_d_abs





# individual changes
frq(data_DE$position_change1) # 42.09 % changed by at least 1 SP
frq(data_DE$position_change2) # 17.42 % changed by at least 2 SP


# direction of opinion change
data_DE <- data_DE %>%
  mutate(diff_pre_post = position_post_cont - position_pre_cont,
         opinion_formation = ifelse(
           abs_position_change == 0, "Stability", ifelse(
             diff_pre_post > 0 & position_pre_cont>=6, "Polarization", ifelse(
               diff_pre_post < 0 & position_pre_cont>=6, "Depolarization", ifelse(
                 diff_pre_post > 0 & position_pre_cont <= 5, "Depolarization", ifelse(
                   diff_pre_post < 0 & position_pre_cont <= 5, "Polarization", NA))))),
         opinion_formation = as.factor(opinion_formation))

frq(data_DE$opinion_formation) # 25.05 % depolarized; 17.04 polarized; 57.91 stability



# 3. Randomization check --------------------------------------------------

# check wether random split into 5 treatment groups worked


frq(data_DE$group_treatment) # cases are equally distributed among 5 groups

# age
ggplot(data_DE, aes(age, group_treatment, group = group_treatment))+
  geom_violin(draw_quantiles = 0.5)+
  geom_jitter(alpha = 0.1)
describeBy(data_DE$age, group = data_DE$group_treatment)
# visual: no differences

# paired t-test
TukeyHSD(aov(age ~ as.factor(group_treatment), data = data_DE))
# no differences between 5 groups regarding age


# education
sjt.xtab(data_DE$education_rec, data_DE$group_treatment, show.col.prc = T)
# no differences between 5 groups regarding education


# pre survey position
ggplot(data_DE, aes(position_pre_cont, group_treatment, group = group_treatment))+
  geom_violin(draw_quantiles = 0.5)+
  geom_jitter(alpha = 0.1)
describeBy(data_DE$position_pre_cont, group = data_DE$group_treatment)
# visual: mean in group 4 (information only) slightly lower than in other groups

# paired t-test
TukeyHSD(aov(position_pre_cont ~ as.factor(group_treatment), data = data_DE))
# differences between 4 (information) and 1 (contestatory)
# differences between 4 (information) and 1 (open communication)
# group 4 (information only) in pre survey slightly more in favor of protecting health


# knowledge
ggplot(data_DE, aes(index_corona_knowledge, group_treatment, group = group_treatment))+
  geom_violin(draw_quantiles = 0.5)+
  geom_jitter(alpha = 0.1)
describeBy(data_DE$index_corona_knowledge, group = data_DE$group_treatment)
# visual: no differences

# paired t-test
TukeyHSD(aov(index_corona_knowledge ~ as.factor(group_treatment), data = data_DE))
# no differences between 5 groups regarding knowledge


# leftright positioning
ggplot(data_DE, aes(leftright, group_treatment, group = group_treatment))+
  geom_violin(draw_quantiles = 0.5)+
  geom_jitter(alpha = 0.1)
describeBy(data_DE$leftright, group = data_DE$group_treatment)
# visual: no differences

# paired t-test
TukeyHSD(aov(leftright ~ as.factor(group_treatment), data = data_DE))
# no differences between 5 groups regarding leftright positioning


# 4. Regression models ----------------------------------------------------


# Study 1 Model A ---------------------------------------------------------


# set opinion stability as reference category
levels(data_DE$opinion_formation)
data_DE$opinion_formation <- relevel(data_DE$opinion_formation, ref = "Stability")

# levels of treatment group
data_DE$group_treatment <- as.factor(data_DE$group_treatment)
levels(data_DE$group_treatment) <- c("Contestatory", "Collaborative", "Open communication", "Information only", "No information (placebo)")
frq(data_DE$group_treatment)
# reorder factor
data_DE$group_treatment <- factor(data_DE$group_treatment, levels = c("No information (placebo)", "Information only", "Contestatory", "Collaborative", "Open communication"))


# missing values
frq(data_DE$opinion_formation) # 2124 valid opinion change; 8 NA
frq(data_DE$group_treatment) # 2132 valid group_treatment; 0 NA
frq(data_DE$strong_prior) # 2125  valid strong_prior; 7 NA
frq(data_DE$index_corona_knowledge) # 2132 valid opinion change; 0 NA
frq(data_DE$leftright) # 2128  valid leftright; 4 NA
frq(data_DE$trust_gov) # 2116  valid trust_gov; 16 NA
frq(data_DE$trust_exp) # 2127  valid trust_exp; 5 NA
frq(data_DE$index_cognition) # 2131  valid index_cognition; 1 NA
frq(data_DE$index_evaluation) # 2130  valid index_evaluation; 2 NA
frq(data_DE$index_accuracy) # 2132  valid index_accuracy; 0 NA


# descriptives continuous independent variables
describe(data_DE$index_corona_knowledge)
describe(data_DE$leftright)
describe(data_DE$trust_gov)
describe(data_DE$trust_exp)
describe(data_DE$index_cognition)
describe(data_DE$index_evaluation)
describe(data_DE$index_accuracy)

# listwise deletion
data_DE <- data_DE %>% 
  filter(!is.na(opinion_formation) & !is.na(group_treatment) & !is.na(strong_prior) &
           !is.na(index_corona_knowledge) & !is.na(leftright) & !is.na(trust_gov) &
           !is.na(trust_exp) & !is.na(index_cognition) & !is.na(index_evaluation) &
           !is.na(index_accuracy))

nrow(data_DE) # 2100 complete cases


# estimate model
de_pol <- multinom(opinion_formation ~ group_treatment + strong_prior + index_corona_knowledge +
                     leftright + trust_gov + trust_exp + index_cognition +
                     index_evaluation + index_accuracy, data = data_DE) 

de_pol_N <- nrow(de_pol$residuals) 
de_pol_N # 2100 complete cases

# coef plot study 1 model A
coef_plot_study1_mA <- plot_model(de_pol, type = "est", show.values = T, transform=NULL,
                                  auto.label = F, vline.color = "black", 
                                  value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                  value.offset =.35, order.terms = c(1:24),
                                  axis.labels= c("Accuracy motivation", "Neeed for Evaluation",
                                                 "Need for Cognition", "Trust experts",
                                                 "Trust government", "Leftright",
                                                 "Corona knowledge", "Strong prior",
                                                 "Open Communication", "Collaborative",
                                                 "Contestatory","Information Only"))+
  labs(title = "(A) All Participants", caption = paste0("N = ", de_pol_N))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))
coef_plot_study1_mA[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study1_mA[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_de_pol <- ggemmeans(de_pol, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_de_pol$x <- factor(pred_de_pol$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study1_mA <- pred_de_pol %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text = element_text(size =7))


study1_plotA <- coef_plot_study1_mA/prob_plot_study1_mA





# Study 1 Model B ---------------------------------------------------------


# filter dataset, keep only pre survey pro health respondents
frq(data_DE$pos_pre_group) # N = 1424

data_DE_health <- data_DE %>% 
  filter(pos_pre_group == 1)

nrow(data_DE_health) # 1424 pro health complete cases


# estimate model
de_pol_health <- multinom(opinion_formation ~ group_treatment + strong_prior + index_corona_knowledge +
                            leftright + trust_gov + trust_exp + index_cognition +
                            index_evaluation + index_accuracy, data = data_DE_health) 

de_pol_health_N <- nrow(de_pol_health$residuals) 
de_pol_health_N # 1424 complete cases

# coef plot study 1 model B
coef_plot_study1_mB <- plot_model(de_pol_health, type = "est", show.values = T, transform = NULL,
                                  auto.label = F, vline.color = "black", 
                                  value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                  value.offset =.35, order.terms = c(1:24),
                                  axis.labels= c("Accuracy motivation", "Neeed for Evaluation",
                                                 "Need for Cognition", "Trust experts",
                                                 "Trust government", "Leftright",
                                                 "Corona knowledge", "Strong prior",
                                                 "Open Communication", "Collaborative",
                                                 "Contestatory","Information Only"))+
  labs(title = "(B) Pro Health", caption = paste0("N = ", de_pol_health_N))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))+
  theme(axis.text.y = element_blank())
coef_plot_study1_mB[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study1_mB[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_de_pol_health <- ggemmeans(de_pol_health, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_de_pol_health$x <- factor(pred_de_pol_health$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study1_mB <- pred_de_pol_health %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text.y = element_blank())


study1_plotB <- coef_plot_study1_mB/prob_plot_study1_mB


# Study 1 Model C ---------------------------------------------------------


# filter dataset, keep only pre survey pro freedom respondents
frq(data_DE$pos_pre_group) # N = 676

data_DE_free <- data_DE %>% 
  filter(pos_pre_group == 0)

nrow(data_DE_free) # 676 pro freedom complete cases


# estimate model
de_pol_free <- multinom(opinion_formation ~ group_treatment + strong_prior + index_corona_knowledge +
                          leftright + trust_gov + trust_exp + index_cognition +
                          index_evaluation + index_accuracy, data = data_DE_free) 

de_pol_free_N <- nrow(de_pol_free$residuals) 
de_pol_free_N # 676 complete cases

# coef plot study 1 model C
coef_plot_study1_mC <- plot_model(de_pol_free, type = "est", show.values = T, transform = NULL,
                                  auto.label = F, vline.color = "black", 
                                  value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                  value.offset =.35, order.terms = c(1:24),
                                  axis.labels= c("Accuracy motivation", "Neeed for Evaluation",
                                                 "Need for Cognition", "Trust experts",
                                                 "Trust government", "Leftright",
                                                 "Corona knowledge", "Strong prior",
                                                 "Open Communication", "Collaborative",
                                                 "Contestatory","Information Only"))+
  labs(title = "(C) Pro Civil-Liberties", caption = paste0("N = ", de_pol_free_N))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))+
  theme(axis.text.y = element_blank())
coef_plot_study1_mC[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study1_mC[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_de_pol_free <- ggemmeans(de_pol_free, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_de_pol_free$x <- factor(pred_de_pol_free$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study1_mC <- pred_de_pol_free %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text.y = element_blank())


study1_plotC <- coef_plot_study1_mC/prob_plot_study1_mC



# show Fig complete
(coef_plot_study1_mA + coef_plot_study1_mB + coef_plot_study1_mC) /
  (prob_plot_study1_mA + prob_plot_study1_mB + prob_plot_study1_mC)
ggsave("study1_regression_opinion change_logit.png", height = 7, width = 12, dpi = 300)





# 5. Model Appendix -------------------------------------------------------


# estimate model A
de_pol_append <- multinom(opinion_formation ~ group_treatment, data = data_DE) 

de_pol_append_N <- nrow(de_pol_append$residuals) 
de_pol_append_N # 2100 complete cases

# coef plot study 1 model A
coef_plot_study1_mA_a <- plot_model(de_pol_append, type = "est", show.values = T, transform= NULL,
                                    auto.label = F, vline.color = "black", 
                                    value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                    value.offset =.35, colors = "#377EB8",
                                    axis.labels= c("Open Communication", "Collaborative",
                                                   "Contestatory","Information Only"))+
  labs(title = "(A) All Participants", caption = paste0("N = ", de_pol_append_N))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))
coef_plot_study1_mA_a[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study1_mA_a[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_de_pol_append <- ggemmeans(de_pol_append, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_de_pol_append$x <- factor(pred_de_pol_append$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study1_mAa <- pred_de_pol_append %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text = element_text(size =7))


study1_plotA_append <- coef_plot_study1_mA_a / prob_plot_study1_mAa


# estimate model B
de_pol_health_append <- multinom(opinion_formation ~ group_treatment, data = data_DE_health) 

de_pol_health_N_append <- nrow(de_pol_health_append$residuals) 
de_pol_health_N_append # 1424 complete cases

# coef plot study 1 model B
coef_plot_study1_mB_a <- plot_model(de_pol_health_append, type = "est", show.values = T, transform = NULL,
                                    auto.label = F, vline.color = "black", 
                                    value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                    value.offset =.35, col = "#377EB8",
                                    axis.labels= c("Open Communication", "Collaborative",
                                                   "Contestatory","Information Only"))+
  labs(title = "(B) Pro Health", caption = paste0("N = ", de_pol_health_N_append))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))+
  theme(axis.text.y = element_blank())
coef_plot_study1_mB_a[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study1_mB_a[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_de_pol_health_append <- ggemmeans(de_pol_health_append, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_de_pol_health_append$x <- factor(pred_de_pol_health_append$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study1_mBa <- pred_de_pol_health_append %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text.y = element_blank())


study1_plotB_append <- coef_plot_study1_mB_a / prob_plot_study1_mBa



# estimate model C
de_pol_free_append <- multinom(opinion_formation ~ group_treatment, data = data_DE_free) 

de_pol_free_N_append <- nrow(de_pol_free_append$residuals) 
de_pol_free_N_append # 676 complete cases

# coef plot study 1 model C
coef_plot_study1_mC_a <- plot_model(de_pol_free_append, type = "est", show.values = T, transform = NULL,
                                    auto.label = F, vline.color = "black", 
                                    value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                    value.offset =.35, col = "#377EB8",
                                    axis.labels= c("Open Communication", "Collaborative",
                                                   "Contestatory","Information Only"))+
  labs(title = "(C) Pro Civil-Liberties", caption = paste0("N = ", de_pol_free_N_append))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))+
  theme(axis.text.y = element_blank())
coef_plot_study1_mC_a[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study1_mC_a[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_de_pol_free_append <- ggemmeans(de_pol_free_append, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_de_pol_free_append$x <- factor(pred_de_pol_free_append$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study1_mCa <- pred_de_pol_free_append %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text.y = element_blank())


study1_plotC_append <- coef_plot_study1_mC_a / prob_plot_study1_mCa



# Fig. appendix complete
(coef_plot_study1_mA_a + coef_plot_study1_mB_a + coef_plot_study1_mC_a) /
  (prob_plot_study1_mAa + prob_plot_study1_mBa + prob_plot_study1_mCa)
ggsave("study1_regression_opinion change_appendix_logit.png", height = 7, width = 12, dpi = 300)



# descriptives study 1
frq(data_DE$strong_prior)
describe(data_DE$index_corona_knowledge)
describe(data_DE$leftright)
describe(data_DE$trust_gov)
describe(data_DE$trust_exp)
describe(data_DE$index_cognition)
describe(data_DE$index_evaluation)
describe(data_DE$index_accuracy)


# 1) Bivariate relationship between communication formats and justification rationality
data_DE_treat <- data_DE %>% filter(data_DE$group_treatment %in% c("Contestatory", "Collaborative", "Open communication"))
data_DE_treat$group_treatment <- droplevels(data_DE_treat$group_treatment)
tab_xtab(data_DE_treat$rj_codes, data_DE_treat$group_treatment, show.col.prc = T, show.summary = F, var.labels = c("Justification rationality", "Treatment"))


# 2) Bivariate relationship between communication formats and Constructive Politics
tab_xtab(data_DE_treat$proposal_final, data_DE_treat$group_treatment, show.col.prc = T, show.summary = F, var.labels = c("Constructive politics", "Treatment"))


# 3) Bivariate relationship between justification rationality and constructive politics
tab_xtab(data_DE$rj_codes, data_DE$proposal_final, show.col.prc = T, show.summary = F, var.labels = c("Justification rationality", "Constructive politics"))


# 4) Bivariate relationship between justification rationality and opinion change
tab_xtab(data_DE$opinion_formation, data_DE$rj_codes, show.col.prc = T, show.summary = F, var.labels = c("Opinion formation", "Justification rationality"))


# 5) Bivariate relationship between constructive politics and opinion change
tab_xtab(data_DE$opinion_formation, data_DE$proposal_final, show.col.prc = T, show.summary = F, var.labels = c("Opinion formation", "Constructive politics"))












# 6. Regression with weights ----------------------------------------------


# Study 1 Model A ---------------------------------------------------------


# estimate model
de_pol <- multinom(opinion_formation ~ group_treatment + strong_prior + index_corona_knowledge +
                     leftright + trust_gov + trust_exp + index_cognition +
                     index_evaluation + index_accuracy, data = data_DE, weights = data_DE$wght) 

de_pol_N <- nrow(de_pol$residuals) 
de_pol_N # 2100 complete cases

# coef plot study 1 model A
coef_plot_study1_mA <- plot_model(de_pol, type = "est", show.values = T, transform=NULL,
                                  auto.label = F, vline.color = "black", 
                                  value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                  value.offset =.35, order.terms = c(1:24),
                                  axis.labels= c("Accuracy motivation", "Neeed for Evaluation",
                                                 "Need for Cognition", "Trust experts",
                                                 "Trust government", "Leftright",
                                                 "Corona knowledge", "Strong prior",
                                                 "Open Communication", "Collaborative",
                                                 "Contestatory","Information Only"))+
  labs(title = "(A) All Participants", caption = paste0("N = ", de_pol_N))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))
coef_plot_study1_mA[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study1_mA[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_de_pol <- ggemmeans(de_pol, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_de_pol$x <- factor(pred_de_pol$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study1_mA <- pred_de_pol %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text = element_text(size =7))


study1_plotA <- coef_plot_study1_mA / prob_plot_study1_mA





# Study 1 Model B ---------------------------------------------------------



# estimate model
de_pol_health <- multinom(opinion_formation ~ group_treatment + strong_prior + index_corona_knowledge +
                            leftright + trust_gov + trust_exp + index_cognition +
                            index_evaluation + index_accuracy, data = data_DE_health, weights = de_pol_health$wght) 

de_pol_health_N <- nrow(de_pol_health$residuals) 
de_pol_health_N # 1424 complete cases

# coef plot study 1 model B
coef_plot_study1_mB <- plot_model(de_pol_health, type = "est", show.values = T, transform = NULL,
                                  auto.label = F, vline.color = "black", 
                                  value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                  value.offset =.35, order.terms = c(1:24),
                                  axis.labels= c("Accuracy motivation", "Neeed for Evaluation",
                                                 "Need for Cognition", "Trust experts",
                                                 "Trust government", "Leftright",
                                                 "Corona knowledge", "Strong prior",
                                                 "Open Communication", "Collaborative",
                                                 "Contestatory","Information Only"))+
  labs(title = "(B) Pro Health", caption = paste0("N = ", de_pol_health_N))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))+
  theme(axis.text.y = element_blank())
coef_plot_study1_mB[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study1_mB[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_de_pol_health <- ggemmeans(de_pol_health, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_de_pol_health$x <- factor(pred_de_pol_health$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study1_mB <- pred_de_pol_health %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text.y = element_blank())


study1_plotB <- coef_plot_study1_mB / prob_plot_study1_mB


# Study 1 Model C ---------------------------------------------------------



# estimate model
de_pol_free <- multinom(opinion_formation ~ group_treatment + strong_prior + index_corona_knowledge +
                          leftright + trust_gov + trust_exp + index_cognition +
                          index_evaluation + index_accuracy, data = data_DE_free, weights = data_DE_free$wght) 

de_pol_free_N <- nrow(de_pol_free$residuals) 
de_pol_free_N # 676 complete cases

# coef plot study 1 model C
coef_plot_study1_mC <- plot_model(de_pol_free, type = "est", show.values = T, transform = NULL,
                                  auto.label = F, vline.color = "black", 
                                  value.size = 2.5, dot.size = 1.5, line.size = .5, 
                                  value.offset =.35, order.terms = c(1:24),
                                  axis.labels= c("Accuracy motivation", "Neeed for Evaluation",
                                                 "Need for Cognition", "Trust experts",
                                                 "Trust government", "Leftright",
                                                 "Corona knowledge", "Strong prior",
                                                 "Open Communication", "Collaborative",
                                                 "Contestatory","Information Only"))+
  labs(title = "(C) Pro Civil-Liberties", caption = paste0("N = ", de_pol_free_N))+
  theme_sjplot()+
  scale_y_continuous(limits = c(-3,3))+
  theme(axis.text.y = element_blank())
coef_plot_study1_mC[["layers"]][[1]][["aes_params"]]$size <- 0.3
coef_plot_study1_mC[["layers"]][[1]][["aes_params"]]$alpha <- 0.5

# predicted probabilities
pred_de_pol_free <- ggemmeans(de_pol_free, terms = "group_treatment", typical = "mean", condition = c(strong_prior =0, gender_male = 0))
pred_de_pol_free$x <- factor(pred_de_pol_free$x, levels = c("Open communication", "Collaborative", "Contestatory", "Information only", "No information (placebo)"))

prob_plot_study1_mC <- pred_de_pol_free %>%
  ggplot(aes(predicted, x))+
  geom_point()+
  geom_segment(aes(x=conf.low, y =x, xend=conf.high, yend=x))+
  facet_grid(~response.level, scales = "free")+
  xlab("")+
  ylab("")+
  theme_bw()+
  scale_x_continuous(n.breaks = 4, labels = label_percent(accuracy = 1))+
  theme(axis.text.y = element_blank())


study1_plotC <- coef_plot_study1_mC / prob_plot_study1_mC


# show Fig complete
(coef_plot_study1_mA + coef_plot_study1_mB + coef_plot_study1_mC) /
  (prob_plot_study1_mA + prob_plot_study1_mB + prob_plot_study1_mC)
ggsave("study1_regression_opinion change_logit_wght.png", height = 7, width = 12, dpi = 300)

